Mass-radius relation of Newtonian self-gravitating Bose-Einstein condensates 
with short-range interactions: II. Numerical results 
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We develop the suggestion that dark matter could be a Bose-Einstein condensate. We determine 
the mass-radius relation of a Newtonian self-gravitating Bose-Einstein condensate with short-range 
interactions described by the Gross-Pitaevskii-Poisson system. We numerically solve the equation of 
hydrostatic equilibrium describing the balance between the gravitational attraction and the pressure 
due to quantum effects (Heisenberg's uncertainty principle) and short-range interactions (scattering). 
We connect the non- interacting limit to the Thomas- Fermi limit. We also consider the case of 
attractive self-interaction. We compare the exact mass-radius relation obtained numerically with 
the approximate analytical relation obtained with a Gaussian ansatz. An overall good agreement is 
found. 



I. INTRODUCTION 

Several recent astrophysical observations of distant 
type la supernovae have revealed that the content of 
the universe is made of about 70% of dark energy, 25% 
of dark matter and 5% of baryonic (visible) matter [l| . 
Thus, the overwhelming preponderance of matter and 
energy in the universe is believed to be dark i.e. unob- 
servable by telescopes. The dark energy is responsible 
for the accelerated expansion of the universe. Its origin 
is mysterious and presumably related to the cosmological 
constant. Dark energy is usually interpreted as a vacuum 
energy and it behaves like a fluid with negative pressure. 
Dark matter also is mysterious. The suggestion that dark 
matter may constitute a large part of the universe was 
raised by Zwicky in 1933. Using the virial theorem to 
infer the average mass of galaxies within the Coma clus- 
ter, he obtained a value much larger than the mass of 
luminous material. He realized therefore that some mass 
was "missing" in order to account for observations. This 
missing mass problem was confirmed later by accurate 
measurements of rotation curves of disc galaxies 0] ■ 
The rotation curves of neutral hydrogen clouds in spi- 
ral galaxies measured from the Doppler effect are found 
to be roughly flat (instead of Keplerian) with a typical 
rotational velocity Vrx> 200km/s up to the maximum 
observed radius of about 50 kpc. This mass profile is 
much more extended than the distribution of starlight 
which typically converges within ^ 10 kpc. This im- 
plies that galaxies are surrounded by an extended halo 



of dark matter whose mass M{r) 



^jG increases lin- 



early with radius [56| . This can be conveniently modeled 
by an isothermal self-gravitating gas the density of which 
scales asymptotically as 

The nature of dark matter (DM) is one of the most 
important puzzles in modern physics and cosmology. A 
wide "zoology" of exotic particles that could form dark 
matter has been proposed. In particular, many grand 
unified theories in particle physics predict the existence 
of various exotic bosons (e.g. ultra-light bosons, ax- 
ions, scalar neutrinos, neutralinos) that should be present 
in considerable abundance in the universe and comprise 



(part of) the cosmological missing mass 0,1. Although 
the bosonic particles have never been detected in accel- 
erator experiments, they are considered as leading can- 
didates of dark matter and might play a significant role 
in the evolution and in the structure of the universe. 

If dark matter is made of bosons, they probably have 
collapsed through some sort of Jeans instability to form 
compact gravitating objects such as boson stars. Boson 
stars were introduced by Kaup 0] and Ruffini & Bonaz- 
zola [l3| in the 1960s. Early works on boson stars [Tl| - [l3 | 
were motivated by the axion field that was proposed as 
a possible solution to the strong CP problem in QCD. 
However, these works consider mini boson stars, like ax- 
ion black holes, with unrealistic small masses. These mini 
boson stars could play a role, however, if they exist in the 
universe in abundance or if the axion mass is extraordi- 
nary small [5^ leading to macroscopic objects with a 
mass Mxaup comparable to the mass of the sun (or even 
larger) [T^. For example, axionic boson stars could ac- 
count for the mass of MACHOs (between 0.3 and 0.8 Mq) 
if the axions have a mass m ^ 10~^°eV/c^ ^M- Alter- 
natively, the possibility that dark matter could be made 
of self-interacting bosonic particles that have formed bo- 
son stars with stellar masses was proposed by Colpi et 
al. [itI and their cosmological formation was discussed 
by Bianchi et al. ^ and Madsen & Liddle [ll|. The 
presence of even a small repulsive self-interaction -jAI^I^ 
between bosons can considerably increase the mass of the 
boson stars allowing therefore to consider larger particle 
masses. For example, for m ~ IGeV/c^ and A ~ 1, this 
mass is of the order of the solar mass Mq, like in the 
case of white dwarf and neutron stars, whereas, in the 
absence of interaction (A = 0), Mxaup ~ 1O~^^M0 for 
m ^ IGeV/c^. Therefore, (mini) boson stars could be 
the constituents of dark matter halos. 

On the other hand, some authors (20l - [35| have pro- 
posed that dark matter halos themselves could be in the 
form of gigantic self-gravitating Bose-Einstein conden- 
sates (with or without self- interaction [5I]) described by 
a single wave function ^p{r,t). In the Newtonian limit, 
which is relevant at the galactic scale, the evolution of 
this wave function is governed by the Gross-Pitaevskii- 
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Poisson (GPP) system Using the Madelung ^ 

transformation, the GP equation turns out to be equiv- 
alent to hydrodynamic (Euler) equations involving a 
barotropic isotropic pressure due to short-range interac- 
tions (scattering) and an anisotropic quantum pressure 
arising from the Heisenberg uncertainty principle. At 
large scales, quantum effects arc negligible and one re- 
covers the classical hydrodynamic equations of cold dark 
matter (CDM) models which are remarkably successful in 
explaining the large-scale structure of the universe [s^l • 
At small-scales, gravitational collapse is prevented by the 
repulsive scattering or by the uncertainty principle. This 
may be a way to solve the problems of the GDM model 
such as the cusp problem [3M| and the missing satellite 
problem [39l |. 

The approach developed in our series of papers 
H^l is essentially theoretical and aims at a general 
study of the basic equations of the problem. In previ- 
ous works, two important limits have been considered. 
Ruffini & Bonazzola [l^l studied boson stars without 
self-interaction described by the Schrodingcr- Poisson sys- 
tem (in the Newtonian limit). The equilibrium state re- 
sults from a balance between gravitational attraction and 
quantum pressure arising from the uncertainty principle. 
On the other hand, Bohmer & Harko [12] studied self- 
gravitating BECs with short-range interactions described 
by the Gross-Pitaevskii-Poisson system and considered 
the Thomas-Fermi (TF) limit in which quantum pressure 
is negligible. In that case, the equilibrium state results 
from a balance between gravitational attraction and re- 
pulsive scattering (for positive scattering lengths a > 0). 
The TF approximation is valid for GN'^m^a/fi^ » 1 f40| 
where N is the number of particles. Since there is some 
indetermination on the values of the mass m and scat- 
tering length a of the bosons, it may be conceptually in- 
teresting to treat the general case in detail and connect 
these two asymptotic limits. This program was started 
in Paper I [40[ using an analytical approach based on a 
Gaussian ansatz and a mechanical analogy. For the sake 
of completeness, we also considered the case of attractive 
short-range interactions (i.e. negative scattering lengths 
a < 0) and reported the existence of a maximum mass 
]^max above which no equilibrium state exists. We ar- 
gued that attractive self-interaction could accelerate the 
formation of structures in the universe. In this paper, 
we shall test the validity of these analytical results by 
determining the exact steady state of the GPP system 
numerically. 

The paper is organized as follows. In Sec. |TT1 we reca- 
pitulate the basic equations of the problem and specifi- 
cally mention the non-interacting case and the Thomas- 
Fermi limit. In Sec. IIIIl we explain our numerical proce- 
dure to solve the quantum equation of hydrostatic equi- 
librium which is equivalent to the steady state of the GPP 
system. In Sec. llVi we plot the fundamental curves (ra- 
dius, energy, density,...) as a function of the mass in 
the case where the value of the scattering length is pre- 
scribed. Finally, in Sec. |Vl we plot these curves as a 



function of the scattering length for a fixed value of the 
total mass. We also compare the exact relations obtained 
numerically with the approximate analytical relations ob- 
tained in Paper I and find an overall good agreement. In 
Appendix \^ we make numerical applications in order to 
compare our results with real astrophysical objects and 
determine the validity of our assumptions. 



II. THE GROSS-PITAEVSKII-POISSON 
SYSTEM 

A. The Madelung transformation 

A Newtonian self-gravitating BEG is described by the 
Gross-Pitaevskii-Poisson system [H, : 
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Ai/j + m{^ + gNm\ilj\'^)ilj, (1) 



(2) 



with g = ATrah^/m^ where a is the s-scattering length 
(we allow a to be positive or negative) [i^. We write 
the wave function in the form ij{r,t) = A(r, t)e''^''''*''/'' 
where A and S are real, and make the Madelung (36j 
transformation 



p = Nm\'ijj\'^ NmA^, 



1 



-VS*, 



(3) 



where p(r, t) is the density field and u(r, t) the veloc- 
ity field. The total mass of the configuration is M = 
Nm = J pdr. We note that the flow is irrotational since 
V X u = 0. With this transformation, it can be shown 
that the Gross-Pitaevskii equation ^ is equivalent to 
the barotropic Euler equations with an additional term 
called the quantum potential (or quantum pressure). In- 
deed, one obtains the set of equations 



dp 

dt 



V • (pu) - 0, 



(4) 



^ + (u.V)u = --Vp-V<i>-lvQ, (5) 
at p m 



with 



and 



2m y/p ' 



2Tmh 2 
P = 5-p 



(6) 



(7) 



The pressure relation ([7]) corresponds to a polytropic 
equation of state of the form 



p = Kp'', 7 = 1 



(8) 
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with polytropic index n = 1 (i.e. 7 = 2) and polytropic 
constant K = 2Trah^ /m^. 

Remark: If dark matter is treated as a coUisionlcss fluid 
(as in standard models [13, the evolution of this 

fluid is fundamentally described by the Vlasov-Poisson 
system. The hydrodynamic equations without pressure 
(h = p = 0) can be obtained if we assume the exis- 
tence of a single fluid velocity at every spatial position 
and they cease to be valid after the first time of cross- 
ing where multi-streaming generates a range of parti- 
cles velocities through a given point [i^. In that case, 
the evolution of the Vlasov-Poisson system becomes very 
complex and exhibits phenomena of nonlinear Landau 
damping and phase mixing [l^. The derivation of fluid 
equations in that context becomes difficult (it amounts 
to closing the infinite hierarchy of Jeans equations) even 
if we may argue that a form of local thermodynamical 
equilibrium and an effective pressure arise as a result of 
violent relaxation [131 . Alternatively, if dark matter is 
a self- gravitating BEC, the hydrodynamic equations 
([5]) are exact for all times since they are rigorously equiv- 
alent to the GP equation ([1]). In that case, the pressure 
p is due to short-range interactions and non-ideal effects, 
not to thermal effects [i^ . There is also a quantum pres- 
sure. These pressure terms prevent multi-streaming and 
the formation of singularities (caustics) , and regularize 
the dynamics at small scales [4^. Therefore, BEC dark 
matter behaves like a fluid, contrary to standard CDM 
that is essentially collisionless. 



C. Hydrostatic equilibrium 

The time-independent solution PT|) can also be ob- 
tained from the Euler equations (jH)-® since they are 
equivalent to the GP equation. The steady state of the 
quantum Euler equation ([5]) obtained by taking dt = 
and u = satisfies 



h p , 



2™2 y ^ J 



(13) 



This is similar to the condition of hydrostatic equilib- 
rium with an additional quantum potential. It describes 
the balance between pressure due to short-range interac- 
tions (scattering), quantum pressure (Hcisenberg princi- 
ple) and gravity. This equation is equivalent to Eq. (jlip . 
Indeed; integrating Eq. (fT3| using Eq. ([7]), we obtain Eq. 
PT|) where the eigenenergy E appears as a constant of in- 
tegration. Combining Eq. ([T3)) with the Poisson equation 
([2|), we obtain the fundamental equation of hydrostatic 
equilibrium with quantum effects 

For the polytropic equation of state ([7]), we get 



47rar ^ r ^ f^Vp\ , ^ 



(15) 



B. The time independent GP equation 

If we consider a wave function of the form 

V(r,i) = ^(r)e-'^^, (9) 
we obtain the time-independent GP equation 



2m 



+ m($ + .giVmi/'^)-0 = (10) 



where -(/'(r) = A{y) is real and p{y) — Nmil}'^{r). The 
foregoing equation can be rewritten 



m$ -|- mgp — 



2m yfp 



or, equivalently, 



m$ -I- mgp + Q = E. 



(11) 



(12) 



In paper I, we have shown that this equation can be ob- 
tained from a variational principle, by minimizing the 
energy functional Etot at fixed mass M. Combined with 
the Poisson equation ([2]), we obtain an eigenvalue equa- 
tion for the wave function ip(r) where the eigenvalue E 
is the energy. In the following, we shall be interested in 
the fundamental eigenmode corresponding to the small- 
est value of E. For this mode, the wave function ip{r) is 
spherically symmetric and has no node so that the den- 
sity profile decreases monotonically with the distance. 



There are three important limits to consider. 

The non-interacting case corresponds to 17 = a = 0. 
This is the situation first studied by Ruffini & Bonazzola 
[l0| and revisited by Membrado et al. [i^ with another 
method. In that case, Eq. ^TE\\ reduces to 



2m2 



m 



AirGp. 



(16) 



It describes the balance between attractive gravity and 
repulsive quantum pressure (Heisenberg principle). This 
equation must be solved numerically. It is found [l^, HI] 
that the density profile decays smoothly to infinity and 
that the radius containing 99% of the mass is given by 
i?99 = 9.9h^/GMm'^. 

For a > 0, we can make the Thomas-Fermi approxima- 
tion which amounts to neglecting the quantum potential. 
This is the limit considered by Bohmcr & Harko [32|. In 
that case, Eq. ([15]) becomes 



Gm^ 
Ap+—^p = 0. 
ah 



(17) 



It describes the balance between attractive gravity and 
repulsive short-range interactions (scattering). This 
equation is equivalent to the Lane-Emden equation for 
a polytropc of index n = 1 [^. It has the analyt- 
ical solution p{r) — (pq-R/t) sin (yrr/i?) where R = 



4 



T^iaT? jG-nt'Y^'^ is the radius of the configuration (inde- 
pendent on the mass M) and po = ttM /AR^ is the central 
density. 

For a < 0, the TF approximation leads to collapse since 
it only keeps the effects of attractive gravity and attrac- 
tive short-range interactions (scattering). If we want to 
obtain equilibrium states (for M < M^ax we have 

to solve the complete Eq. ([T5|) expressing the balance 
between repulsive quantum pressure (Heisenberg princi- 
ple) and attractive short-range interactions (scattering) 
and gravity. We could also consider the non- gravitational 
limit. In that case, Eq. ([T5|) reduces to 
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(18) 



It describes the balance between repulsive quantum pres- 
sure (Heisenberg principle) and attractive short-range in- 
teractions (scattering). However, in Paper I, we have 
shown that such equilibria are unstable. 



III. NUMERICAL SOLUTION 

A. The fundamental differential equation 

The structure of a self-gravitating BEC with short- 
range interactions is determined by the condition of hy- 
drostatic equilibrium (|15l) . In the general case, this equa- 
tion has to be solved numerically. To that purpose, it is 
convenient to introduce dimensionless variables [i^. We 
introduce the lengthscale 



b = 



2GMrr? ' 



(19) 



which basically corresponds to the radius of a self- 
gravitating BEC in the absence of short-range interac- 
tion [loj . We then define the dimensionless position x 
and the dimensionless density profile 7i(x) by 



r inb^ 



(20) 



We also introduce the dimensionless eigenenergy e and 
the dimensionless gravitational potential (/)(x) through 
the relations 



-E. (^(x) - ^^^(r). 



Finally, we introduce the dimensionless parameter 
2Na AGAPma 



X = 



(21) 



(22) 



It is equal to the ratio between N times the scattering 
length a and the typical radius & of a self-gravitating 
BEC without interaction. The non-interacting limit cor- 
responds to X <C 1 and the TF limit corresponds to 1 



In terms of these dimensionless variables, the equation 
of hydrostatic equilibrium (jlSp can be rewritten 



A ( I - yAti = n. 



(23) 



The density must satisfy the normalization condition 



J n{x) dx ~ An. 



(24) 



Finally, the normalized eigenenergy e can be obtained 
from the steady state equation pT|) leading to 



(P + Xn ^ = e, 



(25) 



where n(x) is the solution of Eq. ([25]) and (/)(x) is given 
by 



47r / X — x' 



(26) 



equivalent to the Poisson equation ([2]). For a spherically 
symmetric distribution, Eq. (j23|) reduces to the ordinary 
differential equation 

6n'3 3n"'n' 2n"^ 
+ 





lOn'n' 


h -n'" - 




X 


nx 


In'^n" 


3n'4 


•n? 





2nx{n" + ^] -2n^ (27) 



where the prime ' denotes d/dx. The normalization con- 
dition (1^^ takes the form 



+00 



nx dx = 1, 



and the steady state equation (pS]) becomes 



e — (j> + xn- 



2n nx 4m? 



(28) 



(29) 



The density behaves near the origin like n{x) ~ riQ 4- 
n2x'^ A- .... Taking the limit x — in Eq. ()29p . and using 
the field equation (pS)) . we find that 



3fl2 

n{x)x dx + x''^o • (30) 

no 



The procedure to determine the equilibrium state is now 
clear. We first have to solve the differential equation 
(P7)) with the boundary conditions n(0) — uq, n'{0) = 
n"'{0) = and n"(0) = 2ri2. The constants no and 712 
have to be determined so as to yield a physical density 
profile at infinity and satisfy the normalization condition 
(PS]). Finally, the eigenenergy e is given by Eq. (PH]). 

To solve this problem, it is convenient to introduce 
the function f{x) = n{x)/nQ and make the change of 
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variables X 



1/4 



X and ^ 



1/2 



X- We thus have to 



solve the ordinary differential equation 



10/'/" 6/'-^ 3/"'/' 2/" 



/2 f/l 



fx 



J2 J3 

with the boundary conditions 



- 2/M if" 



f 

2/' 

'Y 



I 

= 2/2, (31) 



/(0) = 1, /'(0) = /"'(0)=0, /"(0) = 



2^2 
3/2 



(32) 



For given (jl, the constant A2 is determined so as to yield a 
physical density profile at infinity (see below) . Then, the 
constant ng is determined by the normalization condition 
which becomes 



1 



1/4 



+ 00 



f{X)X^dX. 



(33) 



Once A2 and no are known, we can obtain the coeffi- 



cient 712 = ^2'T-o^^/2, the parameter x — /^/"-o'^ ^n'^ ^^'^ 



'0 

normalized eigenencrgy 



1/2 



f{x)xdx + ^a^J^ 







^. (34) 
no 



we have determined (A2)* with a precision of the or- 
der of 10-12, i.e. (^2)* ^ -0.612386937160, to get 
X„iax — 18.66. The normalized central density and the 
normalized eigenenergy are found to be no ~ 6.911 lO"'^ 
and e ~ —8.13810-2. Finally, the radius x containing 
99% of the mass is Agg ~ 19.89. These values are consis- 
tent with those obtained in 14911. 
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FIG. 1: Illustration of the shooting problem for determining 
the optimal value of (yl2)«. This example corresponds to x = 
/X = but it is representative of more general cases. 
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B. The numerical procedure 

The numerical procedure that was followed to deter- 
mine the proper value of A2, denoted (A2)*, can be un- 
derstood from the example shown in Fig. [TJ If we solve 
the differential equation (|3T|) with A2 > (A2)* (full hues), 
the density profile reaches a minimum fmin at ^ — X^ 
before increasing indefinitely. As A2 approaches (A2)* 
from above, the point Xc is pushed further and fur- 
ther away while fmin decreases. Therefore, the diver- 
gence of the density occurs at larger and larger radii. 
On the other hand, if we solve the differential equation 
with A2 < (A2)* (dashed lines), the density profile de- 
creases until a point X = Xc at which the program 
breaks down because the density achieves too small val- 
ues (< 10-11). As A2 approaches (^2)* from below, the 
point Xc is pushed further and further away so that the 
"break down" occurs at larger and larger radii. Ideally, 
if we could start exactly from (^2)*, the point Xc would 
be rejected to +00 and the density profile would gen- 
tly decrease towards zero at infinity. In practice, it is 
impossible to obtain the "exact" value of (A2)*. Fur- 
thermore, it is shown in Fig. [5] that the value of Xc 
increases very slowly (logarithmically) with the distance 
\A2 — (^2)*! to the exact value, so that a huge precision 
on the value of A2 is needed to obtain large values of 
Xc- Let us call X^ax the largest value of Xc that we 
have been able to obtain numerically. For x — A* — Oi 




FIG. 2: Evolution of Xc as A2 approaches (A2),. We find 
that the convergence is logarithmic Xc oc ln|j42 — (^2)*! so 
that a huge precision on the value of A2 is needed to obtain a 
large value of Xmax ■ This example corresponds to x = fJ- ~ 
but it is representative of more general cases. 



More generally, we have solved this shooting problem 
for several values of fi. The series of equilibria, that 
will be studied in the following sections, is parameter- 
ized by fi taking values between —00 and -I- 00. Figure 
|3] shows that the parameter x is not a monotonic func- 
tion of fi in the region corresponding to negative scat- 
tering lengths (a < 0). Indeed, x starts from 0- when 
fi — > —00, reaches a minimum value x* — —4.100 at 
/i* ~ —1.000, returns to at /.t = and finally increases 
to +00 when fi — > -|-oo. The TF limit corresponds to 
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FIG. 3: Evolution of the parameter x as a function of /i. 



/i — > +00, the non-interacting limit corresponds to /i = 
and the non- gravitational limit corresponds to /x — — oo. 
The non-monotonicity of x(/x) is associated with an insta- 
bility. This is related to the Poincare theorem (see, e.g., 
50l . 5l| ) since the parameter x plays the role of the mass 
for a given scattering length a. Since the stable equilib- 
rium configurations arc minima of energy at fixed mass 
[40| . a direct application of the Poincare theorem (sec 
Sec. IIV Pp implies that a change of stability occurs at a 
turning point of mass, hence of x- Since we know that 
the system is stable in the TF limit (because it is equiva- 
lent to a polytrope of index 7 = 2 larger than the critical 
index 7c = 4/3 [52l|). we conclude from the Poincare the- 
orem that all the configurations with /i > /i* are dynam- 
ically stable (S) while all the configurations with /i < /i* 
arc unstable (U). In particular, the system is unstable in 
the non-gravitational limit. The same conclusions have 
been reached in Paper I based on the analytical model. 



C. The dimensionless parameters 



is fixed. Using the definitions of Paper I, we normal- 
ize the scattering length by aq = /GM^m, the radius 
by Rq = h^/GMm^, the density by pq = M / R^q = 
G^M'^m^/h^, the total energy by Eq = GAP/Rq = 
G'^lVPm^/n^ and the eigenenergy by E'q = GMm/Rq = 
G^M^m^/h^. In terms of these variables, we obtain 
a/ag = x/4, r/Rg = x/2, R99/RQ = A99/2, E/E'q = 2e 
and p{r)/pQ = (2/T:)n[x). The TF limit is valid for 
a 3> ag and the non- interacting limit for a <^ aq. 

• Let us now consider the case where the scattering 
length a is fixed. We normalize the mass by Ma = 
h/^Gm\al the radius by Ra = {\a\n^ /Grr/Y/"^ , the 
density by Pa ~ Ma/R^ ~ Gra^ / o?h? ^ the total energy 
by Ea = GMl/Ra = h{Gm)^/'^ /\a\^/'^ and the eigenen- 
ergy by E'^ = GMaTn/ Ra = Gm? /\a\. In terms of these 
variables, we have M/Ma ~ -^1x1/2, r/Ra = x/y^jxl, 

Rg^/Ra = Agg/^, E/E'a = e|x|/2 and p{r)lpa = 
X^n(x)/87r. The TF limit is vahd for M > Ma and 
the non- interacting limit for M ^ Ma- 

We can now plot the fundamental curves in the two 
situations described above. The procedure is the follow- 
ing. For a given value of /i, we determine successively 
(A2)*, no, 712, X e. Then, we obtain the density pro- 
file n(x). From the density profile n{x), we obtain the 
value A99 of the radius x containing 99% of the mass. Fi- 
nally, by varying p, between — c» and -|-oo, we obtain the 
fundamental curves (parametrized by p) for a fixed value 
of the scattering length or for a fixed value of the total 
mass by using the parameters defined above. The quan- 
tities that appear in these curves (mass, radius, density, 
energy,...) are dimensionless so that we are not required 
to specify the characteristics of the bosons (that are nec- 
essarily uncertain). In this sense, our approach is very 
general and can model different systems such as (mini) 
boson stars and galactic halos. In Appendix]^ we make 
numerical applications in order to compare our results 
with real astrophysical objects. 



As we have seen, the structure of the problem depends 
on a single control parameter x given by Eq. (j22p . There 
are, however, two ways to present the results. In the first 
case, we assume that the scattering length a is fixed and 
we study how the physical parameters like the radius, 
the energy, the density,... depend on the mass. This 
is certainly the most relevant representation for astro- 
physical problems. In particular, we shall determine the 
mass-radius relation M{R) of boson stars as was done 
in the past for white dwarfs @ and neutron stars [53j . 
In the second case, we assume that the mass M is fixed 
and we study how the physical parameters depend on the 
scattering length. In particular, we shall determine the 
radius versus scattering length relation R(a). This rep- 
resentation can also be of interest. Therefore, we shall 
treat the two situations successively. In each case, the 
physical parameters (mass, radius, scattering length,...) 
must be properly normalized as explained below. 

• Let us first consider the case where the mass M 



IV. THE FUNDAMENTAL CURVES FOR A 
FIXED VALUE OF THE SCATTERING LENGTH 



A. The mass-radius relation 

The mass-radius relation for a fixed value of the scat- 
tering length is plotted in Figs. HI O [5] and [7] for pos- 
itive and negative scattering lengths (Rgg corresponds 
to the radius containing 99% of the mass). In the non- 
interacting case a = 0, or for M — > when a ^ 0, we 
find 



Rgg = 9.946 



GMm^' 



in agreement with previous works [lol |4E 
limit valid for M — )■ +00 (when a > 0) [3^ 



(35) 

. In the TF 

we have the 



7 



analytical result 



R 



99 



2.998 



ah 



1/2 



(36) 



For a > 0, there exists an equilibrium state for all the 
values of the mass M and the configurations are stable 
(S), see Figs. |l]and[5j The radius ([55)1 represents the 
minimum radius achievable by the system [40]. For a < 0, 
there exists a maximum mass 



1.012- 



\/|a|Gm 
corresponding to the radius 



R 



99 



(37) 



(38) 



There is no equilibrium state with Af > Mmax, see Figs. 
[S] and [71 In that case, the system is expected to col- 
lapse and form a black hole [S^]. For M < M„ 



the 

right branch (i? > i?*) is stable (S) and the left branch 
{R < R^) is unstable (U). As indicated previously, the 
change of stability occurs at the turning point of mass in 
agreement with the Poincare theorem (see Paper I and 
Sec. HVDl) . 

The mass-radius relation M{R) for a fixed value of 
the scattering length a has been studied in Paper I by 
using a Gaussian ansatz. This leads to the approximate 
analytical relation 
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FIG. 4: M as a function of 7?99 for fixed a > 0. The mass 
is normalized by Ma and the radius by Ra- In the non- 
interacting limit M — >■ 0, we get M ~ 9.946/-R99. In the 
TF limit M -foo, we obtain _R99 2.998. The system is 
always stable. The dotted line corresponds to the approxi- 
mate analytical mass-radius relation M = 2aR/{yE? — GttQ 
(with i?99 — 2.38 1 67i?) based on the Gaussian ansatz [4Q]. 
The radius is given as a function of the mass by i? = 
{(j/vM){l 4- a/I -I- (StiC^vA'P / cj'^). In the non-interacting limit 
M 0, we get M ~ 2a/uR i.e. M^"-'^'" ~ 8.955/7?99 
and in the TF limit M -l-oo, we get R {<6-nC,/vY''^ i.e. 
^Gaufls _^ ^ -^25. The analytical mass-radius relation has the 
same qualitative shape as the numerical curve and provides 
a good quantitative agreement in the non-interacting limit. 
The agreement is less good in the TF limit where the density 
profile sensibly differs from a Gaussian. 



M 



2cr 



Gm'^R 



-1 67T{^ah-^ 



(39) 



with cr = 3/4, C = l/(27r)3/2, u = 1/^2^. We can also 
express the radius as a function of the mass as 



V GMra^ 



1± Wl 



(40) 



with + when a > and ± when a < 0. This relation is 
compared with the exact mass-radius relation in Figs. [H 
and [Hi for positive and negative scattering lengths. The 
analytical model gives the same scalings as Eqs. ([55|. 
([M| . ([571) and ^ with the prefactors 8.955, 4.125, 1.085 
and 4.125 respectively. We see that the agreement be- 
tween the analytical relation and the numerical one is 
qualitatively correct in all cases. It is also quantitatively 
good, except in the TF limit. This is because, in the TF 
limit, the density has a compact support that is poorly 
represented by a Gaussian distribution. 



B. The mass-energy relation 
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FIG. 5: M as a function of ii99 in log-log plot for fixed a > 0. 
The dashed line corresponds to M ~ 9.946/i?99. 



lengths. In the non-interacting case a = 0, or for A/ — > 
when a 7^ 0, we obtain 



The eigenenergy E is plotted as a function of the mass 
M in Figs. [51 and [HI for positive and negative scattering 



E = -0.1628 



G'^APw? 



(41) 
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FIG. 6: M as a function of Rgg for fixed a < 0. The mass 
is normalized by Ma and tlie radius by Ra- In the non- 
interacting limit M — ^ and R99 — >■ +00, we get M ~ 
9.946/J?99. There exists a maximum mass Mmax — 1.012 
corresponding to a radius Rqq ~ 5.5. The configurations with 
small radius 7? < i?* (i.e. fj, < /j., where fit = —1.000 cor- 
responds to the mass peak) are unstable. The dotted line 
corresponds to the approximate analytical mass-radius rela- 
tion M = 2aR/{uR^ + QtvC) (with ^99 = 2.381677?) based 
on the Gaussian ansatz pj| . The radius is given as a func- 
tion of the mass by J? = (cr/i/M)(l ± ^1 - 6n(uA'P/cr^). In 
the non-interacting limit Rgg — >■ +cxj, we get M ~ 2(j/i>R 
i.e. Af^'""'" ~ 8.955/i?99 and in the non-gravitational limit 
R ^ Q (unstable), we get M ~ cri?/(37rC) i. e. A/ G''"''^ ~ 
0.5262i?99. On the other hand, M^^T" = (^/V^t^ ^ 1-085 
and i??"""" = (67rC/;/)i/2 i.e. (ij^g)''"""" ~ 4.125. The ana- 
lytical mass-radius relation has the same qualitative shape as 
the numerical curve and provides a good quantitative agree- 
ment for any values of the mass. 



For a < 0, the eigenenergy corresponding to the point of 
maximum mass is 

= -0.36-^. (43) 




M 



FIG. 8: iJ as a function of M for fixed a > 0. The eigenen- 
ergy is normalized by jS„ and the mass by Ma- In the non- 
interacting limit M 0, we get E ~ -0.1628M^ and in 
the TP limit M -s- +00, we obtain E ~ -0.3183M. These 
asymptotes are represented in dashed lines. We have also rep- 
resented in dotted line the analytic expression obtained from 
the Gaussian ansatz. It is given hy E ^ a/R'^ + AttCM/R^ - 
2vM/ R where R is related to M by the equation given in the 
caption of Fig. 21 In the non-interacting limit M — >■ 0, we 

get E {■iv'^/4,a)M'^ yielding jc;'^''"'" 0.1592A/2 and in 

the TF limit M +00, we get E ~ -(4i/^''V(3(67rC)^''^))A^ 

yielding E^""""" ' 0.3071A/. The agreement is very good 

for all values of mass and we hardly see the difference between 
the two curves. 
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FIG. 7: M as a function of i?99 in log-log plot for fixed a < 0. 
The dashed line corresponds to M ~ 9.946/i?99. 



In the TF limit valid for M — > -|-oo (when a > 0), we 
have the analytical result [iOl : 



1 G3/2m5/2Af 
E = 



(42) 



The eigenenergy E can be estimated analytically from 
the Gaussian ansatz by using the results of Paper I (see 
in particular Sec. III.B). Using Eqs. (1-77) and (1-83), 
we get 



E = cj- 



2Trah'^M GMm 
2C — 2i^- 



R 



(44) 



where the mass and the radius are related to each other 
by Eq. (p9)) . We can therefore express E as a function 
of R alone 



E = -cr- 



li?2 



- 1 



(45) 



Then, eliminating the radius between M[R) and E(R), 
we obtain as a function of M in parametric form. Al- 
ternatively, we can substitute the analytical expression 
of R{M) given by Eq. (|40]) in Eg. (1341) . In the non- 
interacting case, using E = 3W/2N [40] and Eqs. (1-83) 
and (1-93), we obtain 



E 



1^ ¥ ■ 



(46) 
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FIG. 9: £ as a function of M for fixed a < 0. The eigenen- 
ergy is normalized by E'^ and the mass by Ma- In the non- 
interacting limit M 0, we get E ~ -0.1628M^. The 
eigenenergy corresponding to the maximum mass A'Irnax — 
1.012 is Et ~ —0.36. We have also represented in dotted line 
the analytic expression obtained from the Gaussian ansatz. 
It is given hy E = a/R^ - 4tv(; M/R^ - 2vM/R where R is 
related to M by the equation given in the caption of Fig. 
|6l In the non-interacting limit M — >■ with i? — >■ 0, we get 

E {2,v^/Aa)M'^ yielding ^;GaKas 0.1592M2 and in the 

non-gravitational limit M — >■ with i5 — s> — oo, we get E ~ 

-(CTV(3(37rC)^))M-2 yielding £;Ga«s« 0.3927/M^ The 

eigenenergy corresponding to the maximum mass Af^^"^" ~ 



1.085 is E': 



-5criy/(187rC) ~ -0.4166. 



In the TF approximation (for a > 0), using E ~ AW/iN 
(ioj and Eqs. (1-83) and (1-94), we find that 



E : 



4^3/2 g3/2^5/2j^,j 



(47) 



In the non-Kravitational limit (for a < 0), using E 
-eg/SiV and Eqs. (1-83) and (1-96), we get 



E ■ 



3(3^C)' M^a^ 



(48) 



Finally, combining Eqs. dH]), (1-99) and (I-lOO), the 
eigenenergy corresponding to the point of maximum mass 
is found to be 



E, 



bau Gm? 
ISttC \a\ 



(49) 



This returns the scalings of Eqs. (gT]), (gl]) and 
with the prefactors -0.1592, -0.3071 and -0.4166 (in 
the non-gravitational case, the prefactor is —0.3927). We 
see in Figs. [8] and [9] that the agreement with the numer- 
ical curves is fairly good, even in the TF limit. This is 
because the value of the energy is less sensitive to the 
detailed form of the density profile. Close to the max- 
imum mass (for a < 0), the agreement is good but not 
excellent. 

We can also estimate the total energy from the Gaus- 
sian ansatz by using the results of Paper I (see in partic- 



FIG. 10: Etot (obtained with the Gaussian ansatz) as a 
function of M for fixed a > 0. The energy is normal- 
ized by Ea and the mass by Ma- This yields Etot = 
aM/R^ + 2-kCM'^/R^ - uM'^/R where R is related to M by 
the equation given in the caption of Fig. [l] In the non- 
interacting limit M — >• 0, we get Etot ~ — (i/^/4cr)M^ yield- 
0.05307M^ and in the TF limit M -l-oo, 
-(2i.3/V(3(67rC)i/'))M2 yielding E^oT^' ~ 
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FIG. 11: Etot (obtained with the Gaussian ansatz) as a func- 
tion of M for fixed a < 0. The energy is normalized by 
Ea and the mass by Ma- This yields Etot ~ crM/R^ — 
2-kCM'^ I R^ -vM'^ I R where R is related to M by the equation 
given in the caption of Fig. [6l In the non-interacting limit 
Af -)■ with Etot -i> 0, we get Etot ~ -(i/^/4cr) Af^ yield- 
ing E^oT°" ~ -0.05307Af^ and in the non-gravitational limit 
Af ^ with Etot +00, we get Etot ~ [a^ /{3{3-K(:f))M-^ 
yielding E^oT^" ~ 0.3927A/"^ The energy correspond- 
ing to the maximum mass M^l'i'"' ~ 1.085 is Sgr^' = 
-a'^v^'^ /{i{Q-KC,f'^) ~ -0.09045. The curve makes a spike 
at Af = Mmax corresponding to jj, — jj.^- For M < Mmax, 
there are two solutions for the same mass but the stable state 
(/i > fi,) corresponds to the state of lowest total energy Etot, 
as expected (note that it corresponds to the state of highest 
eigenenergy E, i.e. highest chemical potential a, in Fig. [9]). 



ular Sec. III.B). Using Eqs. (1-75) and (1-83), we get 



E, 



tot 



r 



c- 



GAP 

/ 

R 



(50) 
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where the mass and the radius are related to each other 
by Eq. ([M)) . We can therefore express Etot as a function 
of R alone 




1 - 



27rCa?t 



1 - 



(51) 



Then, eliminating the radius between M{R) and EtotiR), 
we obtain Etot as a function of A/ in parametric form. 
Alternatively, we can substitute the analytical expression 
of R{M) given by Eq. (gD]) in Eq. ((5pl) . In the non- 
interacting case, using Etot = NE/3 [40|, we obtain 



Etnt — — - 



4ct 



(52) 



In the TP approximation (for a > 0), using Etot ~ NE/2 
lioll. we find that 



E,, 



2^3/2 ^3/2„^3/2^j2 



3(60i/2 



(53) 



In the non-gravitational limit (for a < 0), using Etot = 
-NE we get 



Et. 



(54) 



Finally, combining Eqs. (1-99) and (I-lOO), the total 

energy corresponding to the point of maximum mass is 
found to be 



^tot - - 



h{Gm) 



1/2 



3(67rC)3/2 |a|3/2 



(55) 



The prefactors are -0.05307, -0.1536, 0.3927 and 
—0.09045 respectively. The total energy Etot obtained 
with the Gaussian ansatz is plotted as a function of the 
mass M in Figs. [10] and [11] for positive and negative 
scattering lengths. 



C. The density profiles 

The central density po vs mass M is plotted in Figs. 
[T2l and [13] for positive and negative scattering lengths. 
In the non-interacting case a = 0, or for M ^> when 
a ^ 0, we find that 



Pq = 4.400 10" 



(56) 



In the TF limit valid for AI 
we have the analytical result 



Po 



+00 (when a > 0) [32], 
1 G3/2Mm9/2 



47r2 a3/2^3 



(57) 



When a < 0, the central density at the point of maximum 
mass is 



(po)* = 0.04 



(58) 
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FIG. 12: po as a function of M for fixed a > 0. The cen- 
tral density is normalized by pa and the mass by Ma- In the 
non-interacting limit M — >• 0, we get po ~ 4.400 lO'^M* and 
in the TF limit M +oo, we get po ~ 2.533 lO'^M. The 
system is always stable. We have also represented in dot- 
ted line the analytic expression obtained from the Gaussian 
ansatz. It is given by po = M/n^^^R^ where R is related 
to A4 by the equation given in the caption of F ig. |4] yield- 
ing po = (;/^Af''/7r^/V^)(l + ^/l + GTcCuA'P/a^)-^. In the 
non-interacting limit M 0, we get po ~ (z^^/8cr^7r^/^)M* 
yielding pff """" ~ 0.003378M'' and in the TF limit M +oo, 
we get po = (u/Gn'^Cf/'^M yielding p^f"""" ~ 0.03456M. 



It is interesting to compare these results with those 
obtained from the Gaussian ansatz. In that case, the 
central density is related to the mass and to the radius 

by 



M 



Po 



7r3/2i?3 ■ 



(59) 



Combining this relation with Eq. (|39|) . we obtain the 
central density as a function of the radius 



2a 



Po 



-3/2 



(60) 



Eliminating the radius between M{R) and po{R), we ob- 
tain the approximate central density versus mass relation 
in parametric form. Alternatively, we can substitute the 
analytical expression of R{M) given by Eq. (j^O]) in Eq. 
(|59p . In the non-interacting case, we obtain 



Po 



8<t3^3/2 ■ 
In the Thomas- Fermi limit (when a > 0), we find that 



(61) 



In the non-gravitational limit (when a < 0), we get 

,3 



(62) 



Po 



37r3/2cy Af2|a|s 



(63) 
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FIG. 13: po as a function of M for fixed a < 0. Tiie cen- 
tral density is normalized by pa and the mass by Ma- In 
the non-interacting limit M — >■ with po — >■ 0, we get po ~ 
4.400 10~^A'/*. The central density at the point of maximum 
mass Mmax — 1.012 is (po)* — 0.04. The configurations with 
high central density po > (po)* (corresponding to p < p*), 
located after the mass peak, are unstable. We have also rep- 
resented in dotted line the analytic expression obtained from 
the Gaussian ansatz. It is given by Po = M/vr^/^i?^ where R 
is related to M by the equation given in the caption of Fig. 
[6] yielding po = '(;/=*M*/7r^/V=')(l ± ^1 - QttCuAP /a^)-^ . 
In the non-interacting limit M — s- with po — >■ 0, we get 
po ~ (j/V8o-^7r3/2)M'' yielding pg^"""' ~ 0.003378M*. In 
the non-gravitational limit M — ^ with po — >■ -foo, we get 
po ~ {a/STv^/'^CfM-'^ yielding pg^'"""' ~ 0.3535/M^ At the 
point of maximum mass, the central density is (po)?"""" = 
aiy/TT^^'^{6TvCf = 0.03751 
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FIG. 15: Same as Fig. [14] for p = 2.86,4.1,5.3,7.3,9.66. In 
the TF limit p — >■ +oo, the normalized density profile tends 
to the asymptotic distribution p(r) = (l/47r'^) sin(r)/r. 




r 



FIG. 16: Density profiles p(r) along the series of equi- 
libria for fixed a < 0. The density is normalized 
by pa and the radius by Ra- The central density in- 
creases as p decreases. The non-interacting limit cor- 
responds to p — >■ and the non-gravitational limit to 
p — ^ — oo. The different curves correspond to p = 
-0.053, -0.55, -1.6, -2.86, -4.4, -7.83, -14.2, -23. The dis- 
tribution becomes unstable when p < p* ~ —1.000 corre- 
sponding to the turning point of mass (see Sec. IIII B|) . 



FIG. 14: Density profiles p(r) along the series of equilibria 
for fixed a > 0. The density is normalized by pa and the 
radius by Ra. The central density increases as p increases. 
The non-interacting limit corresponds to p — and the TF 
limit to p — >■ -l-oo. The different curves correspond to p = 
0.078,0.34,0.69,1. 



Finally, the central density at the point of maximum mass 
is found to be 



(Po)* 



7r3/2(6^C)' a^n^' 



(64) 



This returns the scalings of Eqs. ([5^ . ([57]) and ((55)) 
with the prefactors 0.003378, 0.03456 and 0.03751 (the 
prefactor in the non-gravitational limit is 0.3535). The 
analytical mass-central density relation is compared with 
the exact relation in Figs. [T^] and [T31 The agreement is 
fairly good, except in the TF limit (when a > 0). Close 
to the maximum mass (when a < 0), the agreement is 
good but not perfect. Some density profiles p{r) are rep- 
resented in Figs. [TH [ini [m for positive and negative 
scattering lengths and different values of M. 
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D. Stability analysis and Poincare theorem 

It is interesting to develop an analogy with thermo- 
dynamics. We have seen in Paper I that stable steady 
states of the GPP system, or equivalently of the quantum 
barotropic Euler-Poisson system, correspond to (local) 
minima of the energy functional (see Eqs. (I-56)-(I-64)): 



U2 , 

P — dv -\ 

^2 2m2 



j (VVp)'dr 



1 



p dr+ - / /9$ dr 



(65) 



at fixed mass M. We thus have to study the variational 
problem 



inm{Etot [p, u] 

p,u 



M[p] = M}. 



(66) 



The critical points, cancelling the first order variations of 
constrained total energy, arc given by SEtot ~ aSM — 
where a is a Lagrange multiplier associated with the con- 
servation of mass that can be interpreted as a chemical 
potential. These first order variations lead to the steady 
state equation (jll[) provided that we make the identifi- 
cation a = E/m. This shows that the eigenenergy E can 
be interpreted as a chemical potential. As a result. Figs. 
[5] and [HI give the chemical potential a = dEtot / dM (con- 
jugate quantity) as a function of the mass M (conserved 
quantity). Similarly, in thermodynamics, the statistical 
equilibrium state of a system is obtained by maximizing 
the entropy S at fixed energy E (see, e.g. jsli])- The 
first order variations are given by 5S — PSE = where 
/3 is a Lagrange multiplier associated with the conserva- 
tion of energy that represents the inverse temperature. 
The caloric curve P{E) gives the inverse temperature 
/3 = dS/dE (conjugate quantity) as a function of en- 
ergy E (conserved quantity). Now, using the Poincare 
theory of linear series of equilibria, we know that when 
we plot the conjugate quantity as a function of the con- 
served quantity, a change of stability can only occur at a 
turning point of the conserved quantity or at a bifurca- 
tion point. In the present case, when a < 0, the change 
of stability occurs at the turning point of mass (see Fig. 
[5]). Similarly, in the thermodynamics of self-gravitating 
systems, the change of stability occurs at the turning 
point of energy in the microcanonical ensemble or at the 
turning point of temperature (equivalent to the turning 
point of mass) in the canonical ensemble [5l|. Finally, 
since 5Etot = SM = at the turning point of mass, we 
conclude that the curve Etot{M) presents a cusp at that 
point (see Fig. Ilip . Similarly, in thermodynamics, since 
5S = SE = at the turning point of energy, the curve 
S{E) presents a cusp at that point (see Fig. 4 in [s^). 

Let us recapitulate: For a > 0, the series of equilibria 
containing all the critical points of the variational prob- 
lem (j66p is parametrized by p going from fj, — >■ -|-oo (TF 
limit) to /i — >■ (non- interacting limit). The series of 
equilibria R{M), E{M) or Etot{M) is monotonic. Since 



the system is stable in the TF limit (see Sec. IIIIB[) . we 
conclude from the Poincare theorem that the whole se- 
ries of equilibria is stable for a > 0. By continuity, the 
non-interacting EEC (/i = a = 0) is also stable. For 
a < 0, the series of equilibria is parametrized by p going 
from /i — > (non-interacting limit) to /i — !■ — oo (non- 
gravitational limit). The series of equilibria R{M), E{M) 
or Etot{M) is non- monotonic. There is a turning point of 
mass at M = Mmax corresponding to = ~ —1.000. 
For M > M,nax there is no solution to the variational 
problem (pS)) and for M < Mmax there are two solutions 
with the same mass. Since the system is stable in the 
non-interacting limit p = a = Q (as we have just seen), 
we conclude from the Poincare theorem that the system 
is stable for p > p^, and unstable for p < p^ (this cor- 
responds to high central densities). This leads to the 
stability/instability regions shown in Figsl4lfT3l 



V. THE FUNDAMENTAL CURVES FOR A 
FIXED VALUE OF THE MASS 

A. The radius vs scattering length relation 

The radius vs scattering length relation for a fixed 
value of the total mass is plotted in Figs. [TTlandfTSl f-Rgo 
corresponds to the radius containing 99% of the mass). 
In the non-interacting case a = 0, the radius is given by 
Eq. psp and in the TF limit a — >■ +oo, it is given by 
Eq. ((36)) . For a < 0, there exists a minimum scattering 
length 



= -1.025 



corresponding to the radius 



Rn 



5.6 



GMm^ ' 



(67) 



(68) 



There is no equilibrium state when a < amin- In that 
case, the system is expected to collapse and form a black 
hole. For a > amin, the upper line (i? > i?*) is stable (S) 
and the lower line {R < R^) is unstable (U). 

The radius vs scattering length relation R(a) for a fixed 
value of the mass M has been studied in Paper I by 
using the Gaussian ansatz. This leads to the approximate 
analytical relation PO]) or inversely 



uGMR^ -2a— R. 



(69) 



The analytical model gives the same scalings as Eqs. 
([551). ([HS]). (|57)) and (UHl) with the prefactors 8.955, 4.125, 
1.178 and 4.477 respectively. 



B. The energy vs scattering length relation 

The eigenenergy E vs scattering length relation for a 
fixed value of the total mass is plotted in Figs. [TH] and 
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FIG. 17: i?99 as a function of a for a fixed value of the mass 
M. The radius is normalized by Rq and the scattering length 
by ag. Stable solutions exist for a > amin — —1.025 and their 
radius Rgg > Rqq ~ 5.6 is monotonically increasing with a. 
In the non-interacting case a = 0, -R99 ~ 9.946 and in the TF 
limit a -s- +00, R99 ~ 2.998a^/^. The non-gravitational limit 
corresponds to a — ^ and 7i99 — > but these solutions are 
unstable. The dotted line corresponds to the approximate 
analytical radius vs scattering length relation a — (uR^ — 
2aR) /QttC (with 7?99 = 2.38 1 677?) based on the Gaussian 
ansatz [4y| . The radius can be expressed in terms of the scat- 
tering length as i? = {a/v){l ± + 67r("i'a/cr2). jjj ^jjg 
interacting case a = 0, we get i? ~ 2cr/i^ i.e. R^^^"'' ~ 8.955 
and in the TF limit a +00, we get R ~ (67rC/i^)^''^a^/^ i.e. 
-RgsT""" ~ 4.125a^/^. In the non-gravitational limit R ^ 
we get R ~ (37rC/cr)|a| i.e. i^g"""" ~ 1.900|a|. On the 



other hand a^?""^ 



'1.178 and i?. 



a/v i.e. 



(-Rgg)'^""'"' ~ 4.477. The analytical radius versus scattering 
length relation has the same qualitative shape as the numeri- 
cal curve and provides a good quantitative agreement except 
in the TF limit where the density profile differs sensibly from 
a Gaussian. 




0.01 0.1 



FIG. 18: i?99 as a function of a in log-log plot for a fixed 
value of the mass M. The dashed line corresponds to i?99 ~ 
2.998a^''^ valid in the TF fimit. 



[20I In the non-interacting case a = 0, the eigenenergy is 
given by Eq. (j41[) and in the TF limit a — > +00, it is given 
by Eq. (j42p . For a < 0, the eigenenergy corresponding 



to the minimum scattering length is 



-0.35- 



(70) 
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FIG. 19: i? as a function of a for a fixed value of the mass 
M. The eigenenergy is normalized by Eq and the scatter- 
ing length by aq. In the non-interacting case a = 0, we 
have E ~ —0.1628 and in the TF limit a — >• +00, we get 
E ~ — 0.3183/a^'''^. The eigenenergy at the point of mini- 
mum scattering length a^nin — —1.025 is iJ* ~ —0.35. We 
have also represented in dotted line the analytic expression 
obtained from the Gaussian ansatz. It is given by _B = a/R? + 
^■K(^a/R? — 2v/R where 7? is related to a by the equation given 
in the caption of Fig. [T71 This yields E ^ -{o + 4:i'R)/3R^ 
or, equivalently, E 
+ 



= -(i^73cr)(5 ± 4^/l + 6-K(iya/a^)/{l ± 
In the non-interacting case a = 0, we 
Gauss ^ _o.i592 and in the 
TF limit a +00, we get E ~ -{4v^/'^ /{3{6nCy^^))a~^^'^ 
yielding E^"''^'' ~ -0.3071/a^/^. In the non-gravitational 
limit M with i? ^ 0, we get E ^ (cr3/(3(37rC)^))a"^ 



get E (3;/74cr) yielding 



yielding E 



'0.3927 /a\ The eigenenergy correspond- 



ing to the minimum scattering length a„ 



-1.178 is 



E. 



Gau 



-0.3537. The analytical eigenen- 



ergy versus scattering length relation has the same qualitative 
shape as the numerical curve and provides a good quantitative 
agreement for all values of the scattering length a. 



Using the Gaussian ansatz, the eigenenergy E is given 
by Eq. where the radius and the scattering lengths 
are related to each other by Eq. We can therefore 

express i? as a function of R alone 

E = -( — +4iyGMmR]—. (71) 
\ m J dn'^ 

Then, eliminating the radius between a{R) and E{R) we 
obtain i5 as a function of a in parametric form. Alterna- 
tively, we can substitute the analytical expression of i?(a) 
given by Eq. (|40p in Eq. pi)) . In the non- interacting 
case, we obtain Eq. (j46|). in the TF limit we obtain Eq. 

and in the non- gravitational limit, we obtain Eq. 
(1481) . At the minimum scattering length, we find that 



E^ 



3(7 



(72) 
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FIG. 20: Eigenenergy E vs scattering length a relation (for 
a > 0) in log-log plot for a fixed value of the mass M. The 
dashed line corresponds to iJ ~ — 0.3183/a^'''^ characterizing 
the TF limit. 
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FIG. 21: Etot (obtained from the Gaussian ansatz) as a func- 
tion of a for a fixed value of the mass M. The total en- 
ergy is normalized by Eg and the scattering length by ag. 
This yields Etot = cr/i?^ + 2nC,a/R^ - p/R where R is re- 
lated to a by the equation given in the caption of Fig. 1171 
This yields Etot = (cr — 2vR)/'iR? or, equivalently, Etot ~ 
-(i/V3cr)(l ± 2\/l + 67rCW'^^)/(l ± + 67rCi^o/o-2)2. In 



the non-interacting case a = 0, we get Eu 



V(4a) 



—0.05305 and in the TF limit a — >• -|-oo, we get Etot ~ 
-{2v^''^ /{i{QTvCf'^))a-^''^ ~ -0.1536/^. In the non- 
gravitational limit a — > 0, we get Etot ~ (o""^/(3(37r(")^))a~^ ~ 
0.3927/a^. The total energy at the point of minimum scat- 
tering length is Etot = -v^ j'irj ~ -0.07074. 



af?M 



(73) 



Then, eliminating the radius between a(i?) and EtotiR) 
we obtain Etot as a function of a in parametric form. 
Alternatively, we can substitute the analytical expression 
of i?(a) given by Eq. (gHl) in Eq. ([50]). In the non- 
interacting case, we obtain Eq. ([5^. in the TF limit we 
obtain Eq. (|53p and in the non-gravitational limit, we 
obtain Eq. (|54p . At the minimum scattering length, we 
find that 



El 



3cr 



(74) 



The prefactor is 0.07074. The total energy Etot is plotted 
as a function of the scattering length a in Fig. [51] for a 
fixed value of the total mass A/. 



C. The density profiles 

The central density po vs scattering length a is plotted 
in Figs. [21] and [231 for a given value of the total mass M. 
In the noninteracting case a = 0, it is given by Eq. (|56p 
and in the TF limit a — > -|-c», it is given by Eq. (j57p . 
The central density at the point of minimum scattering 
length is 



(Po)* 



0.04 



(75) 



With the Gaussian ansatz, the central density is re- 
lated to the mass and to the radius by Eq. ([5^]) . This 
equation gives the approximate central density versus ra- 
dius relation for a fixed mass. Eliminating the radius 
between a{R) and po{R), we obtain po as a function of 
a in parametric form. Alternatively, we can substitute 
the analytical expression of R{a) given by Eq. ([40]) in 
Eq. (|59p . In the non- interacting case a = 0, we obtain 
Eq. (|6ip , in the Thomas- Fermi limit a — > +oo we obtain 
Eq. (|62p and in the non-gravitational limit we obtain 
Eq. ([63]) . Finally, the central density at the point of 
minimum scattering length is found to be 



This returns the scalings of Eqs. ^Q, ([32]) and ([TD]) 
with the prefactors -0.1592, -0.3071 and -0.3537 (in 
the non-gravitational case, the prefactor is —0.3927). We 
see in Fig. [19] that the agreement is fairly good, even in 
the TF limit. 

We can also use the Gaussian ansatz to estimate the 
total energy Etot- It is given by Eq. ([SO]) where the radius 
and the scattering lengths are related to each other by 
Eq. ([M]) . We can therefore express Etot as a function of 



(Po)* 



7]-3/2jj3 



(76) 



This returns the scahngs of Eqs. ([55]), ([57]) and ([75]) 
with the prefactors 0.003378, 0.03456 and 0.02703 (the 
prefactor in the non-gravitational limit is 0.3535). The 
analytical curve is compared with the numerical curve in 
Fig. [22] and the agreement is fairly good. Some density 
profiles p(r) are represented in Figs. [21] [23 [21] for a fixed 
value of the total mass and different values of a. 
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FIG. 22: po as a function of a for a fixed value of the mass Af . 
The density is normalized by pq and the scattering length by 
OQ. In the non-interacting case a = 0, we get po ~ 4.400 10"'' 
and in the TF limit a -i> -foo, we get po ~ 2.533 10"Va^''^- 
The central density at the point of minimum scattering length 
amin — —1.025 is (po)* — 0.04. We have also represented in 
dotted line the analytic expression obtained from the Gaus- 
sian ansatz. It is given by po = I/tt^^'^R^ where R is re- 
lated to a by the equation given in the caption of Fig. 1171 
This yields po = {u'^ /Tv'^/^a-^)/{l ± ^1 + G-KCva/a'f. In the 

yield- 



non- interacting case a = we get po ~ {u'' / 



ing pg'"""' ~ 0.003378, in the TF limit a -5> -|-cxj we get 
po = {v/6Tr''(f^''<^~^^' yielding p^""""" ~ 0.03456/0^/^ and 
in the non-gravitational limit a —> with po — > +oo we get 
po ~ (cr/37r^/2(^)3/|^|3 yielding pg'"""'' ~ 0.3535/|a|^ Finally, 
the central density at the point of minimum scattering length 
is (po)f """" = j-K^I'^a^ = 0.02703. The approximate ana- 
lytical curve has the same qualitative shape as the numerical 
curve and provides a good quantitative agreement for all val- 
ues of the scattering length. The most severe deviation occurs 
close to the minimum scattering length value but it remains 
weak. 
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FIG. 23: Po as a function of a (for a > 0) in log-log plot for 
a fixed value of the mass M. The dashed line corresponds to 
Po ~ 2.533 10~V«^''^ characterizing the TF limit. 



FIG. 24: Density profiles p(r) along the series of equilib- 
ria for a fixed value of the total mass M. The density is 
normalized by pq and the radial distance by Rq. The cen- 
tral density decreases as [i, increases. The non-gravitational 
limit corresponds to p — >■ — oo, the non-interacting case to 
p = and the TF limit to p — > -fco. We have repre- 
sented p = —2.86, —2.55, —1.6, —0.55, —0.053 (negative scat- 
tering lengths). The distribution becomes stable when p > 
p. ~ —1.000 corresponding to the turning point of scattering 
length (see Sec. HITBI) . 
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FIG. 25: Same as Fig. [24]for p = 0,0.078,0.34,0.69, 1 (posi- 
tive scattering lengths). 



VI. CONCLUSION 

In this paper, "we have obtained the exact mass-radius 
relation of self-gravitating BECs with short-range inter- 
actions by numerically solving the equation of hydro- 
static equilibrium taking into account quantum effects. 
We have also compared our results "with the approximate 
analytical relation obtained in Paper I from a Gaussian 
ansatz. We have found that the Gaussian ansatz al"ways 
provides a good qualitative agreement "with the exact so- 
lution, and that the agreement is quantitatively good in 
many cases. This gives us confidence to extend our ap- 
proach to more general situations. This is interesting be- 
cause analytical methods allo"w to explore a "wider range 
of parameters than numerical methods and to obtain a 
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FIG. 26: Same as Fig. [24]for = 4.1,4.4,5.3,7.3 (positive 
scattering lengths). 



more complete picture of the problem in parameter space. 
This will be the object of Paper III. 



Appendix A: Numerical applications 

In this Appendix, wc try to relate our general results 
to real astrophysical systems and test the various approx- 
imations made in our study. To facilitate the comparison 
with astrophysical data, we express the radius _R, the to- 
tal mass M, the scattering length a and the boson mass 
m in terms of Mq, kpc, fm and eV/c^. 



1. Fundamental quantities 

Let us first consider the case where the scattering 
length a is fixed. The characteristic radius Ra = 
{\a\TL^ / G'wP'Y/'^ and the characteristic mass Ma = 

hi 



Gmlal can be written as 



AI„ o/i / Ifm 

— = 1.537 10"^^ — 



1/2 



leV/c2 



m 



3/2 



1/2 /I -sr/ 2\ 1/2 

^IcV/c^^ 



(Al) 



(A2) 



For a > 0, the characteristic radius Ra corresponds to 
the minimum radius Rmin = T^Ra of a Newtonian self- 
gravitating BEC, i.e. the radius of the BEC in the TF 
limit, and the characteristic mass Ma corresponds to the 
typical mass separating the TF regime (M ^ Ma) from 
the non-interacting regime (M <C Ma). For a < 0, 
the characteristic mass Ma corresponds to the maximum 
mass Mmax = 1.012Afa and Ra is the corresponding ra- 
dius i?gg = 5.5i?a- 

Let us now consider the case where the total mass M 
is fixed. The characteristic radius Rq = Ti^ / GMra^ and 



the characteristic scattering length ag = fi^ /GM'^m can 
be written as 



8.544 10-- 
Ikpc M 



m 



I«qI 

Ifm 



= 2.36310" 



-68 



M 



leV/c2 



(A3) 



(A4) 



For a > 0, the characteristic radius Rq corresponds 
to the minimum radius R^™ = 9.946i?Q of a Newto- 
nian self-gravitating BEC, i.e. the radius of the BEC 
in the non-interacting limit, and the characteristic scat- 
tering length aq corresponds to the typical scattering 
length separating the TF regime (a 3> ag) from the non- 
interacting regime (a <C ag). For a < 0, the character- 
istic scattering length aq corresponds to the minimum 
scattering length a,„m = — 1.025aQ and Rq is the corre- 
sponding radius i?gg = 5.6i?Q. 

In order to make the correspondence between BECs 
with short-range interactions described by the Gross- 
Pitaevskii equation and scalar fields with a ■|A|0|'* in- 
teraction described by the Klein-Gordon equation, it is 
useful to introduce the dimensionless parameter [40[: 



a 



(A5) 



where A,, = ti/mc is the Compton wavelength of the 
bosons. It can be written as 



— = 5.06810-^— 

Stt Ifm leV/c2 



(A6) 



Using this expression, we can express the results in terms 
of A and m instead of a and m. In particular, wc find 
that Ra = {\\\h^/87rGm'^c)^/^ and Ma = i87:hc/G\X\)^/^ 
leading to 



= 78.08 
Ikpc V 87r \ m 



M^ 

Mr. 



1.09410" 




(A7) 



(A8) 



Finally, we have seen in Paper I that the timcscalc of 
the BEC oscillations (for stable systems with a > 0) or 
the collapse of the BEC (for unstable systems with a < 0) 
corresponds to the dynamical time ~ [R^ /GM)^/^ . It 
can be written 



^ = 1.488 10^^ ( 

Is \lkpc 



3/2 

M 



(A9) 



On the other hand, the maximum circular velocity is 
{vc)max ^ {GM/R)^/"^ = R/to- It can be written 



iVc)r 



1 km/s 



2.07410" 



/Ikpc 



1/2 



M 
M^ 



1/2 



(AlO) 
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2. Orders of magnitude 



3. Validity of the TF approximation 



A specificity of our approacli is to remain as general 
as possible (although restricting ourselves to the Newto- 
nian approximation) so that we have expressed our re- 
sults for an arbitrary value of the scattering length a 
(positive, zero, or negative) and of the boson mass to. 
Therefore, our results can describe different objects, at 
different scales, such as (mini)-boson stars or galactic ha- 
los. In this subsection, and in the following subsections, 
we consider the case of galactic halos and make some nu- 
merical applications in order to check the validity of our 
approximations. To that purpose, we take a radius of or- 
der R = 10 kpc and a total mass of order M = 3 lO^^M© 
js^ . The corresponding dynamical time to is of the order 
of 27 Million years and the maximum circular velocity of 
order 360km/s. 

Let us first assume that the galactic halo can be mod- 
eled by a non-relativistic gas of self-gravitating BECs 
with short-range interactions in the TF approximation. 
In that case, its radius R = nRa is determined by Eqs. 
(|A1[) and (|A7|) . and it is independent on the total mass. 
We have 



R 
Ikpc 



0.01746 



R 



Ikpc 



Vlfm 



48.93 



1/2 



leV/c2 



3/2 



(All) 



(A12) 



We see that the values of a and to given by Bohmer 
and Harko [s^j, namely (a, to) = (1 fm, 14 meV/c^) and 
{a,m) = (10^ fm, 1.44 eV/c^) give the correct order of 
magnitude of the radius. For these values, A/Stt = 
7.095 10"" and A/Stt = 7.297 lO'^ respectively. In fact, 
it is important to realize that the radius i? of a sclf- 
couplcd Newtonian BEC directly determines the ratio 
a/m? or A/to"*. Using Eqs. (|XTT|) and (|M2l) with 
R = 10 kpc, we obtain m^/a = 3.049 10-*^(cV/c2)3/fm 
and m^jX = 23.94 (eV/c^)^ in agreement with the esti- 
mate of Arbey et al. [3l| (they find a larger numerical 
coefficient ^ 50 because they take a smaller halo radius). 

Let us now assume that the galactic halo can be mod- 
eled by a non-relativistic gas of sclf-gravitating BECs 
without short-range interaction. In that case, its radius 
i?99 = 9.946i?Q is determined by Eq. (jA3p leading to 



R 



99 



Ikpc 



84.9810"^^ 



Mr. 



M \ 



/leV/c 



(A13) 



We sec that the mass of the bosons must be very small, of 
the order of to = lO^^^cV/c^, to reproduce the correct 
values of the radius and mass of the cluster. This boson 
mass corresponds to the estimate of Baldeschi et al. [l^l 
and others [2l|, [H, [H-d^ • 

Therefore, a self-interaction can increase the required 
value of the boson mass from to — 10~^^eV/c^ to 
171 = 1 eV / which may be more realistic from a par- 
ticle physics point of view. 



For a > 0, the validity of the TF approximation is 
determined by the dimensionless parameter (40| : 



X 



(A14) 



The TF approximation is valid for x ^ 1 while the non- 
interacting limit corresponds to x ^ 1. We can write 
X = {M/Maf = ajaq. For a fixed value of a, the TF 
approximation is valid for M 3> Ma and for a fixed value 
of Af , the TF approximation is valid for a 2> oq. Equa- 
tion (|JT4)) can be rewritten 



X = 4.23210 



67 



M 



a TO 



MqJ lfmleV/c2' 



or, equivalently. 



X = 8.361 10 



75 



M 
M^ 



8tt' 



(A15) 



(A16) 



Let us make a numerical application. If we want to model 
a galactic halo, then M ~ 3 IO^Mq. Therefore, the TF 
approximation is valid if 



a TO 



Ifm leV/c2 
or, equivalently, 

A 



> 2.62610" 



— > 1.33110 



-99 



(A17) 



(A18) 



This relation clearly shows that the limit A — !■ is dif- 
ferent from the non-interacting case A = 0. Indeed, the 
TF approximation is valid even for a (very) small value 
of A fulfilling the condition (jAl8[) . By contrast, if A = 
strictly, we are in the non-interacting case. In these two 
extreme limits, the physics of the problem is very differ- 
ent (see Appendix lA 21 and Sees. II. E and II. F of Paper 
I). For the values of a and to given by Bohmer and Harko 
j32| . see Appendix I A 21 the condition (|A18p is fulfilled by 
more than 90 orders of magnitude (!) so that the TF 
approximation is perfect. 



4. Validity of the Newtonian approximation 

The Newtonian approximation is valid if the ra- 
dius R of the configuration is much larger than the 
Schwarzschild radius Rs = 2GM/c^, or, equivalently, if 



Rc^ 
G 



For fixed a ^ 0, this condition can be rewritten 

\a\c^ \X\hc 



M R 



with K 



Gm SirGm^ ' 



(A19) 



(A20) 
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The relativity parameter k can be written as 

47 \a\ leV/c2 



K = 7.55410 



Ifm 



m 



or, equivalently, 



(A21) 



Mr, 



= 0.1011 VAfi^ 
Y m 



i?^ 9.272^^ km. 



(A27) 



(A28) 



. ^ 1.491 10- P^l^^/^' 



(A22) 



We also note that k = (| A|/87r)(Afp/TO)^ where Mp = 
ihc/Gf/'^ is the Planck mass. To determine the region in 
the {M,R) plane of Figs. |4] and [6] where the Newtonian 
approximation is valid, it suffices to draw a straight line 
with slope K. The Newtonian approximation is valid be- 
low this line. In general, k is large (see Eqs. (jA21[) and 
IA22[) ) so that the Newtonian approximation is a good 
approximation as soon as a 7^ 0, even if it is (very) small. 
Let us try to be more specific by considering particular 
limits. 

In the TF limit, the radius of the classical self- 
gravitating BEC is of the order Ra- Using Eq. (jA19[) . we 
find that the Newtonian approximation is valid if 



(Gm)3/2 



A 1 



3/2 



(A23) 



Using Eqs. (|XT|) . (|X6)) and (|A20|) . this condition can be 
rewritten 



M 



< 1.16110 
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or, equivalently, 
M 



_Q_\i/2 / leV/c^ 
Ifm/ \ m 



3/2 



< 1.631 10'" 1/ — 




(A24) 



(A25) 



If we want to model a galactic halo, recalling the numer- 
ical values of M, a and m given in Appendix lA 21 we 
see that the term in the left hand side is of order 10^^ 
while the term in the right hand side is of order 10^^. 
Therefore, the condition (|A24p - (jA25l) is fulfilled by 6 or- 
ders of magnitude so that the Newtonian approximation 
is very good. We note that the mass appearing in the 
right hand side of Eq. (|A23[) represents the maximum 
mass of a relativistic self-gravitating BEC with short- 
range interactions in the TF limit (up to a numerical 
factor that can be obtained by solving the general rel- 
ativistic equation of hydrostatic equilibrium (ssf). As 
shown in Appendix B3 of Paper I, it corresponds to the 
maximum mass Mmax = 0M2VXMf>/m'^ found by Colpi 
et al. [l7| . The corresponding radius is still given by 
Ra which can be written R = 0.3836-\/A(Afp/m)Ac or 
R = 6.187GM^ax/c^ m. We have 



Mq ' Vlfm/ V m 



3/2 



(A26) 



For m ^ IGeV/c^, a ^ Ifm and A 1, the maximum 
mass Mmax ^ Mq is of the order of the solar mass, and 
the corresponding radius R is of the order of the kilo- 
meter, like in the case of white dwarf and neutron stars. 
This could describe boson stars with relevant mass. We 
emphasize that the mass M (or the radius R) of the self- 
coupled relativistic boson stars directly determine the ra- 
tio a/m^ or A/m^. Taking M = 1 Mq, we obtain /a ~ 
1.302 (GeV/c2)3/fm and m^/A = 0.01023 (GeV/c^)^. 

In the non-interacting case (a = 0), the radius of the 
classical self-gravitating BEC is of the order of i?Q. Using 
Eq. (|A19p . we find that the Newtonian approximation is 
valid if 



he 



or 



i?> — . (A29) 



Using Eqs. (|A3|) and (|A20p . this condition can be rewrit- 
ten 



M 
M^ 



< 1.33610" 



-10 



leV/c2 



or 



R 
Ikpc 



07 leV/c^ 
> 6.39510-27 



m 



(A30) 



(A31) 



If we want to model a galactic halo, recalling the numeri- 
cal values of M and m given in Appendix lA 21 we see that 
the term in the left hand side of Eq. (jA30p (resp. Eq. 
(jA31[) ) is of order 10^^ (resp. 10) while the term in the 
right hand side is of order 10^^ (resp. 10"^). Therefore, 
the conditions (|A30p and (|A3ip are fulfilled by 3 orders 
of magnitude so that the Newtonian approximation is 
very good. We note that the mass and the radius ap- 
pearing in the right hand side of Eq. (jA29p represent the 
maximum mass and the minimum radius of a relativis- 
tic self-gravitating BEC without short-range interaction. 
As shown in Appendix B.2 of Paper I, it corresponds to 
the Kaup mass Mxaup = 0.633Afp/m and Kaup radius 
i?95 6.03Ac 01 R = 9.526GM^ax/c^- We have 



'^'^''-^ 0.8457 10-1° 



Afp 



i?95 = 14.07^%^ km. 

Mq 



(A32) 



(A33) 



For m ~ IGeV/c^, the Kaup mass MKaup ~ lO-^^A^g 
is very small and irrelevant. The Kaup mass becomes 
of the order of the solar mass if the bosons have a mass 
m ^ lO-i^eV/c^ (leading to a radius of the order of the 
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km). For example, axionic boson stars could account for 
the mass of MACHOs (between 0.3 and 0.8 Mq) if the 
axions have such a small mass [T6j . 

Finally, in the case of attractive short-range interac- 
tions (a < 0), the Newtonian approximation is valid 



if K > {M,nax/Ma)/{R*Qg/Ra) = 0.184, a condition 
that is easily fulfilled in practice (see Eqs. (|A2ip and 
(jA22[) ). As shown in Appendix B.4 of Paper I, we have 
M„,, = 5.073Mp/^\ and R^g = l.l^\{Mp/m)X,. 
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] The maximum mass (|37p that we obtain in the case of 
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